{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Histogram matching algoritm with numpy\n",
    "\n",
    "Pure numpy algorithm to adjust values of a **source** ndarray such that it's histogram \n",
    "matches that of a **reference** ndarray. The source and reference arrays can be any dimension\n",
    "and may differ in shape. They should have the same dtype however."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "source2D = np.array([[10, 20, 40],\n",
    "                   [20, 10, 60],\n",
    "                   [40, 90, 10]])\n",
    "\n",
    "reference2D = np.array([[ 25, 35, 63],\n",
    "                      [63, 25, 90],\n",
    "                      [63, 99, 10]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the original shape for later to reconstruct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 3)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "orig_shape = source2D.shape\n",
    "orig_shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the rest of the analysis, unravel into a 1D array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[10 20 40 20 10 60 40 90 10]\n",
      "[25 35 63 63 25 90 63 99 10]\n"
     ]
    }
   ],
   "source": [
    "source = source2D.ravel()\n",
    "print(source)\n",
    "reference = reference2D.ravel()\n",
    "print(reference)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the set of unique pixel values and their corresponding indices and counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[10 20 40 60 90] [0 1 2 1 0 3 2 4 0] [3 2 2 1 1]\n"
     ]
    }
   ],
   "source": [
    "s_values, s_idx, s_counts = np.unique(\n",
    "    source, return_inverse=True, return_counts=True)\n",
    "\n",
    "print(s_values, s_idx, s_counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at what `s_values` and `s_counts` do. It's simply all unique values in the source and their corresponding counts. So you can see we have three 10s, two 20s, etc. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((10, 3), (20, 2), (40, 2), (60, 1), (90, 1))\n"
     ]
    }
   ],
   "source": [
    "print(tuple(zip(s_values, s_counts)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `s_idx` contains indexes to look up values from `s_values` in order to reconstruct the original array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10, 20, 40],\n",
       "       [20, 10, 60],\n",
       "       [40, 90, 10]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the original source array\n",
    "np.array([s_values[idx] for idx in s_idx]).reshape(orig_shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can do the same with the reference (but we don't need the `r_idx` for this)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[10 25 35 63 90 99] [1 2 1 3 1 1]\n"
     ]
    }
   ],
   "source": [
    "r_values, r_counts = np.unique(reference, return_counts=True)\n",
    "print(r_values, r_counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to calculate the empirical cumulative distribuition, scaled 0 to 1. Each `quantiles` tells us, for each unique `value`, what proportion of the data fall at or below that value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((10, 0.33333333333333331),\n",
       " (20, 0.55555555555555558),\n",
       " (40, 0.77777777777777779),\n",
       " (60, 0.88888888888888884),\n",
       " (90, 1.0))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s_quantiles = np.cumsum(s_counts).astype(np.float64) / source.size\n",
    "tuple(zip(s_values, s_quantiles))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((10, 0.1111111111111111),\n",
       " (25, 0.33333333333333331),\n",
       " (35, 0.44444444444444442),\n",
       " (63, 0.77777777777777779),\n",
       " (90, 0.88888888888888884),\n",
       " (99, 1.0))"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r_quantiles = np.cumsum(r_counts).astype(np.float64) / reference.size\n",
    "tuple(zip(r_values, r_quantiles))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1093625f8>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VFX6x/HPoUrvRVpCUbAgi1IFIRQTbIAodYHAAvJT\nWVfRXZHdxegqCy6iuGIFpO3SRCQISCAwSDEQejH0shCK9BYgZZ7fH3cIISZkSGbmTibP+/XKi8zM\nzbkPk/Dl5NxzzzEiglJKqcCVz+4ClFJKeZcGvVJKBTgNeqWUCnAa9EopFeA06JVSKsBp0CulVIDL\nMuiNMRONMSeNMdtuc8wnxpi9xpgtxpjfebZEpZRSOeFOj/4bICyzF40xTwC1ReQeYDDwhYdqU0op\n5QFZBr2IrAbO3eaQTsBU17HrgFLGmEqeKU8ppVROeWKMvipwJM3jeNdzSiml/IBejFVKqQBXwANt\nxAPV0zyu5nruN4wxurCOUkplg4iY7H6tuz164/rISCTQF8AY0ww4LyInM2tIRPzq4+2337a9htxS\nl9akNeW2upxO4cAB4ZtvhH79hOBgoUIF4fnnhX//W9i2TUhJsf+9OnXlFCGTQ3jqP09x4dqF37ye\nU1n26I0x/wVCgHLGmP8BbwOFrMyWr0RkkTHmSWPMPuAK0D/HVSmlVDaIwP79sHIlOBzWn0lJ0Lq1\n9fHnP8N994HJdt/Y83b8uoNOMzvx/H3PM7LdSPLny+/xc2QZ9CLSy41jhnimHKWUcp8I7Nlza7Ab\nczPY//Y3uPde/wr2tCJ3RzIgcgAfhX1E74d6e+08nhijz9VCQkLsLiFD/liX1uQercl9d1qXCMTF\n3Qz1n36CQoWsUG/XDv7xD6hVK2fB7ov3SkQYtXoU42PH80PPH2haralXz2c8Mf7j9smMEV+eTymV\nuzmdsGOHFeo3gr14cSvYQ0KsP4OD7a7yzlxNusqAyAHsPbuX77t/T9WSWc9GN8YgObgYq0GvlPIb\nKSmwbdvNYF+1CsqUuTXYq1fPshm/FX8xns6zOnNP2XuY2HEiRQoWcevrNOiVUrlWcjJs2XIz2Fev\nhooVb4Z6q1ZQNUBuv1wfv54us7rwcuOXGdZyGOYOxpc06JVSuUZSEmzadDPY16yBatVuXjxt1Qoq\nV7a7Ss+bvm06ry15jYkdJ9Kxbsc7/noNeqWU30pMhA0bbs6K+flnqFnz1mCvUMHuKr0nxZnC8Ojh\nzPllDpE9I3mw4oPZakeDXinlN65fh/Xrbwb7unVwzz03g/2xx6BcObur9I2L1y/Sa24vriRdYU7X\nOZQvWj7bbWnQK6Vsc/WqFeY3pjtu2AD16t28eNqyJZQubXeVvrfv7D46zuhI66DWfPLEJxTMXzBH\n7WnQK6V8JiEB1q69Oca+aRM8+ODNi6ctWkDJknZXaa/lB5fTc25PIlpH8GLjFz3Spga9UsprLl+2\nLpjeCPatW6FBg5vB/uij1rx2Zd0E9VnsZ/zjp38w47kZtKnZxmNta9ArpTzm4kVriuONYN+xAx55\n5OYYe/PmULSo3VX6n8SURF5Z/Aqr/reKyB6R1C5b26Pta9ArpbLt/HnrpqQbwb5rFzRufDPYmzaF\nIu7d05NnnU44zfOzn6dk4ZJM7zKdkoU9P3alQa+UctuZMzeD3eGAffugWbObwd6kCRQubHeVucf2\nk9vpNLMT3R/ozntt3/PKypOgQa+Uuo1Tp6z1YW4E+6FD1rj6jVkxjzxiLQqm7tz8XfMZuGAgH4d9\nzO8f+r1Xz6VBr5RKdfLkrUv2xsdbM2FuXDx9+GEokOfXrM0ZEWHkqpF8vuFz5nWfR+Oqjb1+Tg16\npfKw+Pib4+srV8Kvv1pz128E++9+B/m9M5qQJyUkJTAgcgD7z+7n+x7fU6VEFZ+c1ydBb4zpAHyM\ntfXgRBEZne710sAkoDZwFfiDiPySQTsa9ErlwJEjN3vrK1fCuXPWMgI3xtjr19dg95ajF4/SeWZn\n6pavy4RnJri98qQneD3ojTH5gD1AO+AYEAv0EJFdaY75ALgkIv8wxtQFxotI+wza0qBX6g4cOnRr\nsF++fDPUW7eGBx6AfO7u/KyyLeZoDM/Nfo5XmrzCX1r85Y5WnvSEnAa9O6N1TYC9InLYdcKZQCdg\nV5pj7gf+CSAiu40xwcaYCiJyKruFKZXXpN3v9MZHYqJ/73eaF0zbOo3Xo15nYseJPFP3GbvLyRZ3\ngr4qcCTN46NY4Z/WVqALsMYY0wSoAVQDNOiVykTa/U5vfMCt+53ec48Gu11SnCm8Ff0Wc+PmsiJ8\nBQ9UfMDukrLNU9ffRwHjjDGbgO3AZiAlowMjIiJSPw8JCfHbvSyV8rQb+52mDfa0+52++27O9ztV\nnnHh2gV6fdeLq0lXWT9wPeWK+nbJTYfDgcPh8Fh77ozRNwMiRKSD6/EwQNJfkE33NQeB+iJyOd3z\nOkav8gynE3buvDXYixe/OSMmN+53mhfcWHmybc22fBT2UY5XnvQEX1yMzQ/sxroYexxYD/QUkbg0\nx5QCEkQkyRgzCGghIv0yaEuDXgW0Q4dg/nzrAuqN/U7TBntu3u800CUkJbBwz0KGLB7CuyHvMrjR\nYLtLSuX1i7EikmKMGQJEcXN6ZZwxZrD1snwF3AdMMcY4gZ3AgOwWpFRuc+kSzJkDU6daPfhnn4Vu\n3WD8eKjim2nWKhtEhG0nt7Fk/xKi9kexLn4dDSs3ZPbzs2kd3Nru8jxKb5hSKhtSUmD5cpgyBX74\nweq1h4fDU0/pkgL+7MTlEyzdv5SoA1Es3b+UkoVLElo7lLDaYYQEh1CicAm7S8yQ3hmrlA/FxVk9\n9+nToVIlK9x79AjsfU9zs2vJ11j9v9VE7Y8ian8Uhy8cpl3NdoTWDuXxWo9Ts0xNu0t0iwa9Ul52\n5gzMnGn13o8ehd69rYB/IPfOtgtYIsIvp36xgv1AFGv+t4b6leoTWiuU0NqhNK7amAL5ct9iPxr0\nSnlBUhIsWmT13qOj4cknrXBv104XBfM3pxNOs+zAstRee8H8BQmrHUZo7VDa1mxL6bty/6a1GvRK\neYgIbN5s9dxnzIC6da1w79oVSpWyuzp1Q2JKIj8f+Tm1177nzB5CgkNSe+11ytbx+RIF3qZBr1QO\nHT9ujblPnQpXrkDfvtCnD9T27G5wKptEhL1n96b22FceXkndcnUJrW0Fe/Nqzf1irrs3adArlQ1X\nr1rz3adMgZgY6NLF6r23bKmLhPmDc1fPsfzg8tSpjymSktpjb1+rvc/vVLWbBr1SbhKBNWusnvu3\n31p7o4aHQ+fOuuG13ZKdyaw7ui51OGbnrztpWaNl6tTHeuXrBdxwzJ3QoFcqC4cOWeE+dao1xz08\n3Jo5U7Wq3ZXlbQfOHUgdjllxaAXBpYMJrRVKWJ0wWlRvQeECunntDRr0SmUg/d2qPXpYY++NGumi\nYXa5eP0iKw6uSO21X7p+KbXH3r5WeyoVr2R3iX5Lg14pF71b1b+kOFPYeHxjaq9984nNNK/WPPUi\nav2K9fP0cMyd0KBXeV5cnBXu06dD5cpWz71nT71b1Q5HLhxJ7bEvO7CMKiWqpF5EfSzoMYoW1Ish\n2aFBr/KkjO5W7dsXHnzQ7sryliuJV1h5eCVR+6NYsn8JpxNO83itx1OXGKhaUi+EeIIGvcozbtyt\nOmXKrXertm+vd6v6ilOcbD2xNXXaY+yxWBpVaZTaa294d0PyGZ2f6mka9Cqg6d2q9jt+6XjqcMzS\n/UspW6Rs6jh7SHAIxQsVt7vEgKdBrwKS3q1qn6tJV1n1v1WpF1GPXjxKu1rtUnvtQaWD7C4xz9Gg\nVwFD71a1h4iw49cdqb32tUfW0qBSg9Spj42qNCJ/vvx2l5mn+STojTEdgI+5ucPU6HSvlwOmA3cD\n+YEPRWRyBu1o0Ktb3LhbdcoUmDvXmuceHm7t0qR3q3rPr1d+vWXFx7sK3EVY7TDC6oTRJrgNpe7S\ncTF/4os9Y/MBe7D2jD0GxAI9RGRXmmPeBu4SkbeMMeWx9pitJCLJ6drSoFeA3q3qa9eTr7P2yNrU\n2TEHzh0gJDgkdTnf2mV1TMyfeX3PWKAJsFdEDrtOOBPoBOxKc8wJoL7r8xLAmfQhr9TFi9YaMzfu\nVu3e3brAqnerep6IsPvM7tRgX3V4FfdXuJ/Q2qF88sQnNK3aNOBXfFQ3uRP0VYEjaR4fxQr/tL4G\noo0xx4DiQHfPlKdyu4zuVv3Tn6ypkYV1KROPOnv17C3DMQBhtcPo16Af056dRtkiZW2uUNnFU7OP\n3wK2ikgbY0xtYKkx5iERuZz+wIiIiNTPQ0JCCAkJ8VAJyp9kdLfqRx/p3aqelJSSRMzRmNSLqHGn\n4mgV1IrQ2qG88egb1C1XV5cYyKUcDgcOh8Nj7bkzRt8MiBCRDq7HwwBJe0HWGLMIeF9E1rgeRwNv\nisiGdG3pGH0A07tVvW/f2X2pPXbHIQe1y9ZOXfGxebXmuuJjgPLFxdj8WBdX2wHHgfVATxGJS3PM\nh8BFEXnHGFMJ2AA0EJGz6drSoA8wiYmweLHereotF65dYPnB5am99qtJV29Z8bFCMf0VKS/w5fTK\ncdycXjnKGDMYq2f/lWumzTdADcAA/xSRGRm0o0EfAERg0ybroqrerepZyc5kNhzbkHoRddvJbTxa\n/dHU2TEPVHhAh2PyIL1hSvnMjbtVp0y5ebdq3756t2pOHT5/ODXYlx9cTrWS1VKDvWWNlhQpWMTu\nEpXNNOiV1509Cy+9BEuW6N2qnhR9IJohi4dwJuFM6toxj9d6nLtL3G13acrP+GIevcrD1q611nZ/\n7jnrAmuxYnZXlPuJCB/+/CEf/vwhkzpOIqxOmK74qLxKg15lyOmE0aPh449hwgR45hm7KwoMVxKv\nMCByAPvO7mPdwHXUKFXD7pJUHqBBr37j5Elr7D0hATZsgOrV7a4oMOw7u49nZz1LoyqNWNV/lY69\nK5/R3xfVLZYvh4cfhsaNYcUKDXlPWbx3MS0mteClRi8xqeMkDXnlU9qjVwAkJ8O771rDNFOnWvPg\nVc45xcn7P73PFxu/4Ltu39GiRgu7S1J5kAa9Ij4eevWCggWt+fGVK9tdUWC4eP0ifef15VTCKWIH\nxVKlRBW7S1J5lA7d5HELF8Ijj0BoqDV9UkPeM+JOxdHk6yZUKVGFFeErNOSVrbRHn0clJsJf/wqz\nZsGcOfDYY3ZXFDjmxc1j8A+DGdV+FH9o+Ae7y1FKgz4vOngQevSAihWtoZry5e2uKDCkOFMYsWIE\n07ZNY2GvhTSu2tjukpQCdOgmz5k7F5o2tYI+MlJD3lPOXT3H0zOeZu3RtWx4YYOGvPIr2qPPI65d\ng9dfhx9/tMblG2sOecy2k9voMqsLHet25IPHP6BAPv1npfyL/kTmAXv2QLducO+91lCNrjDpOTN3\nzOSPi//IuA7j6FW/l93lKJUhDfoAN306vPYavPcevPCC7s3qKcnOZN5c+ibf7/6eZX2W0aByA7tL\nUipTGvQB6soV+OMfrUXJli2DBppDHnPqyim6f9udgvkLEjsoVvdiVX5PL8YGoO3brTF4p9Naq0ZD\n3nM2HttIo68b0axaMxb1WqQhr3IFt4LeGNPBGLPLGLPHGPNmBq+/YYzZbIzZZIzZboxJNsaU9ny5\n6nZE4OuvoW1bePNNmDwZihe3u6rAMXnLZDr8pwMfhX3EyHYjyZ8vv90lKeUWd/aMzQfswdoz9hgQ\nC/QQkV2ZHP808KqI/Ga1FN14xHsuXrTG4OPirJug6tWzu6LAkZiSyGs/vsayg8uY130e91e43+6S\nVB6T041H3OnRNwH2ishhEUkCZgKdbnN8T+A3+8Uq79m40VpxskwZiInRkPek45eO03ZKW45eOsr6\nges15FWu5E7QVwWOpHl81PXcbxhjigAdgLk5L01lRQTGjYMnnoB//hM+/xyK6Oq3HrP2yFoaf92Y\nsNphzOs+j1J36bxUlTt5etbNM8BqETmf2QERERGpn4eEhBASEuLhEvKGs2ehf384dszqxdeqZXdF\ngUNE+HLjl4xYMYJvOn3DU/c+ZXdJKo9xOBw4HA6PtefOGH0zIEJEOrgeDwNEREZncOx3wGwRmZlJ\nWzpG7wE39nF9/nmrJ1+okN0VBY5rydd4eeHLrItfx7zu87in3D12l6SUT8boY4E6xpggY0whoAcQ\nmUEhpYDWwPzsFqNuz+m0gv3ZZ+HTT+HDDzXkPenIhSO0+qYVlxIvETMwRkNeBYwsh25EJMUYMwSI\nwvqPYaKIxBljBlsvy1euQzsDS0TkqvfKzbt0H1fvchxy0HNuT4Y2G8obj76B0VuIVQDJcujGoyfT\noZtsWb4c+vSxxuQjIqCA3s/sMSLCxzEfM3rNaKY9O43Haz9ud0lK/UZOh240MvyY7uPqXQlJCQxa\nMIi4U3HEDIwhuHSw3SUp5RUa9H5K93H1rgPnDtBlVhceqvQQa/6whiIFdV6qCly61o0f0n1cvStq\nfxTNJzZnQMMBTOk8RUNeBTzt0fuRtPu4fvsttGxpd0WBRUQYtXoU/17/b+Z0nUOroFZ2l6SUT2jQ\n+4m0+7hu3gzlytldUWC5dP0S/ef35+jFo6wftJ5qJavZXZJSPqNDN34g/T6uGvKetefMHppNbEbZ\nImVZ2W+lhrzKc7RHbyPdx9X7IndHMjByIO+3fZ9BjwyyuxylbKFBbxPdx9W7nOLkHcc7TNoyicie\nkTSr1szukpSyjQa9DXQfV+86f+08vb/rzcXrF9kwaAOVileyuySlbKVB70Np93GNjoaHHrK7osCz\n89eddJ7VmSfqPMGHoR9SMH9Bu0tSynZ6MdZH0u/jqiHveXN2ziFkSggjWo3gkyc+0ZBXykV79F4m\nYi1hMHw4jBkD4eF2VxR4UpwpDI8ezuxfZrOk9xIevvthu0vyW8HBwRw+fNjuMlQmgoKCOHTokMfb\n1aD3orT7uK5apVv8ecOZhDP0mNsDESF2UCzli5a3uyS/dvjwYXRhQf/lrVVTdejGS3QfV+/bfHwz\njb5uRMPKDfmx948a8kplQnv0HiYCn3wC778P48dD1652VxSYpm2dxtCooYx/cjzdHuhmdzlK+TW3\ngt4Y0wH4mJsbj2S0jWAI8BFQEDglIm08WGeuoPu4el9SShJvRL3Bon2LWBG+ggcrPmh3SUr5vSyH\nbowx+YBPgTDgAaCnMaZeumNKAeOBp0XkQSDP9WPXroWGDaFOHVizRkPeG05ePkn7ae3Zd24fsYNi\nNeSVcpM7Y/RNgL0iclhEkoCZQKd0x/QC5opIPICInPZsmf7rxj6uXbpYQzW6j6t3rDu6jkZfNyIk\nKIQFPRdQ+q7SdpekVK7hTtBXBY6keXzU9Vxa9wJljTErjDGxxpg+nirQn508CU88AYsWQWwsPP20\n3RUFpgmbJvDMjGf49IlPeafNO+QzOocgUK1evZoWLVpQunRpypcvz2OPPcbGjRvtLivX89TF2ALA\nw0BboBjwszHmZxHZ56H2/Y7u4+p915Ov88riV/jpfz+xqv8q6pava3dJyosuXbrEM888w5dffknX\nrl1JTExk1apVFC5c2KPnSUlJIX/+/B5t09+50zWKB2qkeVzN9VxaR4ElInJNRM4APwENMmosIiIi\n9cPhcGSjZHslJ8OIEdC7N0yZYq1XoyHvefEX4wmZEsLpq6dZP3C9hnwesGfPHowxdOvWDWMMhQsX\npn379jz44IOICO+99x7BwcFUrlyZfv36cenSJQBWrlxJ9erVb2mrZs2aLF++HIB33nmHrl270qdP\nH0qXLs2UKVNwOp2MHDmSOnXqUKpUKRo3bkx8vBVru3btIjQ0lHLlynHfffcxZ84c374RgMPhuCUr\nc0xEbvsB5Af2AUFAIWALcF+6Y+oBS13HFgW2A/dn0JbkZkePirRqJdK+vcjx43ZXE7h+OvST3D3m\nbhn500hxOp12lxNQ/Pnf4MWLF6V8+fISHh4uixcvlnPnzqW+NnHiRLnnnnvk0KFDcuXKFenSpYv0\n6dNHREQcDodUr179lraCg4MlOjpaREQiIiKkUKFCEhkZKSIi165dkw8++EAeeugh2bt3r4iIbNu2\nTc6ePStXrlyR6tWry5QpU8TpdMqWLVukQoUKEhcX54u3INPvj+v5LPM6sw/3DoIOwG5gLzDM9dxg\n4IU0x7wB7AS2AX/MpB3Pvis+9MMPIpUqibz/vkhKit3VBCan0ymfxHwiFf9VURbvXWx3OQEpq3+D\n1p0gOf/Irl27dkn//v2levXqUqBAAenUqZOcPHlS2rVrJ59//nnqcbt375ZChQpJSkqKW0HfunXr\nW16vW7euLFiw4DfnnzVrlrRq1eqW5wYPHizvvvtu9v9Sd8BbQe/WoIOI/AjUTffcl+kejwHGuPVr\nRC6i+7j6xtWkq/zfwv9jy4kt/DzgZ2qV0fmpdrB7dYS6desyadIkwBrK6d27N6+++irHjx8nKCgo\n9bigoCCSkpI4efKkW+2mH9o5cuQItTKYA3348GFiYmIoW7YsYHWEU1JS6NMnd88v0ekLt3HwIDz2\nGOzaZe3jqiHvHYfPH6blNy1JSkli7R/WasgrAO69917Cw8PZsWMHVapUuWUxtsOHD1OwYEEqVapE\nsWLFSEhISH0tJSWFU6dO3dJW+jVkatSowf79+39zzurVqxMSEsLZs2c5e/Ys586d4+LFi4wfP97D\nfzvf0qDPhO7j6hvRB6JpOqEpvev35j9d/kOxQsXsLknZZPfu3YwdOzb1ouiRI0eYMWMGzZs3p2fP\nnnz00UccOnSIy5cv89e//pUePXqQL18+7r33Xq5du8bixYtJTk7mvffeIzEx8bbnGjBgAH//+9/Z\nt8+aGLh9+3bOnTvH008/zZ49e5g+fTrJyckkJSWxYcMGdu3a5fW/vzfpfJF0dB9X3xARxqwdw9iY\nscx4bgZtaua5FTNUOiVKlGDdunWMHTuWCxcuULp0aZ555hk++OADihUrxrFjx2jVqhXXr1+nQ4cO\nfPLJJwCULFmSzz77jAEDBuB0OvnLX/5CtWq33wB+6NChJCYmEhoaypkzZ6hXrx7z5s2jTJkyREVF\n8dprrzF06FBEhAYNGjB27FhfvAVeY8SHg3LGGPHl+e5U2n1cv/5a93H1lsuJlxkQOYAD5w4wt9tc\napSqkfUXKY8wxugyxX4ss++P6/lsr2GsQzcu06dDixbw4ovWhVcNee/Yd3YfzSc2p1jBYqzqv0pD\nXikfyPNDN7qPq+8s2ruI/vP7E9E6gv9r9H9e22RBKXWrPB30p09Dq1bQpIm1j2vx4nZXFJic4uT9\nn97ny41fMq/7PB6t/qjdJSmVp+TpMfpu3aBaNcjl11n82oVrFwj/PpxTCaf4tuu33F3ibrtLytN0\njN6/6Ri9h82aBdu3w8iRdlcSuOJOxdF0QlOqlKjCivAVGvJK2SRP9uhPnIAGDWDBAmvYRnnevLh5\nDP5hMKPbj6Z/w/52l6NctEfv37zVo89zY/QiMHgwDByoIe8NKc4URqwYwfTt01nYayGNq+qNCErZ\nLc8F/fTp1tIGs2fbXUngOXv1LL3m9uJ6ynViB8VSsVhFu0tSSpHHxujj4627XqdMAQ/vZZDnbT2x\nlcZfN+b+CveztM9SDXnlU59//jmVK1emZMmSnDt3zu5y/E6eGaMXgSefhGbN4O23bSkhYM3YPoNX\nfnyFcR3G0at+L7vLUbfhz2P0wcHB/PrrrxQoUIBixYrRvn17PvvsM0qUKHHbr0tOTqZkyZKsX7+e\nBx/M3RvG66ybHJo4EX79FYYPt7uSwJHsTOb1Ja/ztxV/Y1mfZRryKkeMMSxcuJCLFy+ydetWtm/f\nznvvvZfl1504cYLr169z3333Zeu8TqczW1+Xm+SJoD98GN56yxqyKVjQ7moCw6krpwidFsrOUzuJ\nHRRLg8oZ7hyp1B250ZutWLEiYWFh7Ny5E4DExETeeOMNgoKCuPvuu3nppZe4fv06e/fupV69egCU\nKVOG9u3bA7ffDrB///689NJLPPXUU5QoUQKHw5Fp+3Bzq8KxY8dSqVIlqlatyuTJk1Pbu3btGq+/\n/jrBwcGUKVMmdeE1gJiYGFq0aEGZMmVo2LAhK1eu9Pp7mCF3difB2mFqF7AHeDOD11sD54FNro+/\nZdLOnW234gEpKSJt24r8858+P3XAio2PlRof1ZDhy4ZLckqy3eWoO2DHv0F3pd0V6siRI1K/fv3U\nnZ1effVV6dSpk5w/f14uX74sHTt2lOHDh4uIyKFDhyRfvnyp205mtB1g+fLlU7cD7Nevn5QuXVp+\n/vlnEbG2Frxd+w6HQwoUKCARERGSnJwsixYtkqJFi8r58+dFROSll16SNm3ayPHjx8XpdMrPP/8s\niYmJEh8fL+XKlZMff/xRRESWLVsm5cqVk9OnT2f6HmT2/cHbWwli9fpv7BlbEGvP2HrpjmkNRLrR\nVubfZS8ZP16kaVORpCSfnzogTdo0Scp/UF7m/jLX7lJUNmT1b5AIPPKRHcHBwVKiRAkpUaKEGGOk\nc+fOkuLat7NYsWJy4MCB1GPXrl0rNWvWFBGRgwcPSr58+VKPzWo7wH79+kl4ePgtr9+ufYfDIUWL\nFk1tX0SkYsWKsm7dOnE6nVKkSBHZvn37b/4+o0ePlr59+97yXFhYmEydOjXT98BbQe/O9MomwF4R\nOQxgjJkJdHL18NPyuxWq9u+HESNg9WookOcmknpWYkoir/34GssOLmNlv5XcX+F+u0tSXiBv23uh\ndv78+bRp04affvqJjh07snHjRoKDg0lISOCRRx5JPc7pdKYO86RfHC+z7QD79u2bekzarQVPnTp1\n2/YBypUrR758N0e6ixYtyuXLlzl9+jTXr1/PdFvC2bNns2DBgtQ6kpOTadu2bbbem5xwJ/6qAkfS\nPD6KFf7pNTfGbAHigT+LyC8eqC/bnE7o39+6+OoawlPZdPzScbrO6Uq5ouVYP3A9pe7SNZyVd9wI\n11atWjFkyBDefPNNoqOjKVq0KDt37uTuu7NeRuPGdoBLlizJ9Ji0/zmUL1/+jtpPq3z58tx1113s\n37+f+vURMnZ/AAAU1ElEQVTr/6aOvn378uWXX2by1b7jqYuxG4EaIvI74FPgew+1m23jxllTKv/0\nJ7sryd3WHllL468bE1Y7jHnd52nIK5959dVXWb9+PevXr2fQoEG8+uqrqXvBxsfHExUVlXps2t53\nZtsB7t69O8PzGGOybD8zxhj69+/P0KFDOX78OE6nk5iYGJKSkujduzcLFiwgKioKp9PJtWvXWLly\nJceOHcvJ25It7gR9PJB2d4hqrudSichlEUlwfb4YKGiMKZtRYxEREakfDocje1VnYfdueP99mDwZ\n8uf3yikCnojweezndJ7ZmS+f/pK/t/47+UyemKSlbJJ+CKZ8+fKEh4czevRoRo8eTZ06dWjWrBml\nS5cmNDSUPXv2ZPi1xYsXJyoqipkzZ1KlShWqVKnCsGHDUmfCZCSr9m9X65gxY6hfvz6NGzemXLly\nDBs2DKfTSbVq1Zg/fz4jR46kQoUKBAUFMWbMGLemczocjluyMqeyvGHKGJMf2A20A44D64GeIhKX\n5phKInLS9XkTYLaIBGfQlmR1vpxKSYGWLaF3b3j5Za+eKmBdS77GywtfZl38OuZ1n8c95e6xuyTl\nIf58w5SycVEzEUkxxgwBorB+A5goInHGmMHWy/IV8Lwx5kUgCbgKdM9uQTk1ZgwUKWJtCaju3JEL\nR3hu9nMElw4mZmAMxQvpbixK5XYBtQTCjh3Qpg3ExkJwsNdOE7Achxz0nNuToc2G8sajb+hWfwFI\ne/T+TZcpzkJSEoSHWxuJaMjfGRHh45iPGb1mNNO7TKd9rfZ2l6SU8qCACfpRo6BCBWudeeW+hKQE\nBi0YRNypOGIGxhBcOtjukpRSHhYQQb9lC/z737BpE+hog/sOnDvAs7OepUGlBqz5wxqKFCxid0lK\nKS/I9fPlEhOtIZt//cva6Fu5Z8m+JTSf2JyBDQcypfMUDXmlAliu79G/+y4EBUGau5vVbYgIo1aP\n4t/r/82crnNoFdTK7pKUDwUFBelFdj8WFBTklXZz9ayb2Fh4+mnYuhUqV/ZYswHr0vVL9Jvfj/iL\n8cztNpeqJavaXZJSyg15duORa9esIZtx4zTk3bH79G6aTmhKuSLlWNlvpYa8UnlIrg36ESPg/vuh\nu223ZuUekbsjeeybx3it2Wt89cxXFC6gG+YqlZfkyjH6tWth2jTYtk1n2dyOU5y843iHSVsmEdkz\nkmbVmtldklLKBrku6BMSoF8/GD/emjevMnby8kn6ft+Xq0lXiR0US+XiOr6lVF6V64Zu3noLGjeG\nLl3srsR/RR+I5uGvHqZxlcYsD1+uIa9UHperevQrV8K338L27XZX4p+Sncm8u/JdJmyawNRnp+pS\nBkopIBcF/eXL1o5RX34JZTNc6T5vO3rxKL3m9qJwgcJsGrxJe/FKqVS5Zujmz3+G1q2tefPqVgv3\nLKTRV43oUKcDS3ov0ZBXSt0iV/Toly6FhQt1yCa9xJREhkcPZ/bO2Xzb7Vta1mhpd0lKKT/k90F/\n4QIMGAATJkAp3a401cFzB+kxtwcVi1Vk8+DNlCtazu6SlFJ+yq2hG2NMB2PMLmPMHmPMm7c5rrEx\nJskY47E5MUOHwhNPQGiop1rM/eb+MpemE5rS88GeRPaI1JBXSt1Wlj16Y0w+4FOsPWOPAbHGmPki\nsiuD40YBSzxV3MKFsHy5dWOUsvZyfX3J6/y4/0cW9lpI46qN7S5JKZULuNOjbwLsFZHDIpIEzAQ6\nZXDcH4FvgV89UdjZszB4MEyaBCVKeKLF3G3PmT00m9CMUwmn2PTCJg15pZTb3An6qsCRNI+Pup5L\nZYypAnQWkc8BjyxK8Mor1k1Rbdp4orXcbdrWabSY1IIXG73IrOdnUeouvVihlHKfpy7GfgykHbvP\nNOwjIiJSPw8JCSEkJOQ3x8ybB+vWWTtH5WVXEq8wZPEQYo7GEN03mocqPWR3SUopH3A4HDgcDo+1\nl+V69MaYZkCEiHRwPR4GiIiMTnPMgRufAuWBK8ALIhKZrq0s16M/dQoeesi6A7ZFizv96wSO7Se3\n0/3b7jSp2oRPn/yU4oWK212SUsomOV2P3p2gzw/sxroYexxYD/QUkbhMjv8GWCAi32XwWpZB360b\n1KgBY8a49xcINCLChE0TGL58OB+GfkjfBrp1llJ5XU6DPsuhGxFJMcYMAaKwxvQnikicMWaw9bJ8\nlf5LslvMrFnWTVFTpmS3hdzt4vWLvLDgBeJOx7Gq/yrqla9nd0lKqQDgN1sJnjgBDRrAggXQpInP\nSvIbG45toMe3PXi81uOMDRurm3UrpVJ5fejGkzILehF49llrx6iRI31Wjl8QET5Z9wnvr3qf8U+O\np+sDXe0uSSnlZ7w+dOML06fDgQPW0E1ecvbqWfrP78/xS8eJGRhDrTK17C5JKRWAbF+9Mj4eXn/d\nGpcvnIe2Ml3zvzU0/LIhdcrUYfUfVmvIK6W8xtYevQgMHAgvvwwNG9pZie84xcno1aMZt24cEzpO\n4Ol7dd1lpZR32Rr0EyfCyZMwfLidVfjOjX1cE5IS2PDCBqqVrGZ3SUqpPMC2oZvDh639X6dMgYIF\n7arCd5YfXM7DXz1MkypNWBG+QkNeKeUztvTonU5rjfmhQ6F+fTsq8B3dx1UpZTdbgv6LL+DSJWt7\nwEAWfzGeXt/1olD+QrqPq1LKNj4futm/H0aMsIZsCvjF5E7vWLhnIY989QhhtcN0H1ellK18HrX9\n+1tj8/UC9O7+xJRE/hr9V2btnKX7uCql/ILPg97phFdf9fVZfUP3cVVK+SOfD91Mngz58/v6rN63\naO8i3cdVKeWX/GKtm9zu+13fM/iHwczvMZ9m1ZrZXY5SKsAExFo3udncX+by8qKXWfz7xTx898N2\nl6OUUr+hQZ8Dc3bO4Y+L/8iPvX/kd5V/Z3c5SimVIdsXNcutZu6YySs/vkJUnygNeaWUX3Mr6I0x\nHYwxu4wxe4wxb2bwekdjzFZjzGZjzAZjTFvPl+o//rPtPwxdMpSlfZbqht1KKb/nzp6x+YA9WHvG\nHgNigR4isivNMUVFJMH1eX1gnojUyaCtXH8xdtrWaQyLHsbSPku5v8L9dpejlMoDcnox1p0efRNg\nr4gcFpEkYCbQKe0BN0LepThwOrsF+bPJWybzVvRbRPeN1pBXSuUa7gR9VeBImsdHXc/dwhjT2RgT\nBywCXvFMef5j0uZJ/H3F34nuG62bdiulchWPzboRke+B740xLYFpQN2MjouIiEj9PCQkhJCQEE+V\n4DVfbfyK9356j+V9l3NPuXvsLkcpFeAcDgcOh8Nj7bkzRt8MiBCRDq7HwwARkdG3+Zr9QBMROZPu\n+Vw3Rv957OeMWjOK6L7R1Cn7m8sOSinldb4Yo48F6hhjgowxhYAeQGS6Imqn+fxhgPQhnxuNXz+e\n0WtGsyJ8hYa8UirXynLoRkRSjDFDgCis/xgmikicMWaw9bJ8BTxnjOkLJAJXgO7eLNoXPln3CR/H\nfIyjn4Pg0sF2l6OUUtmma91k4KOfP+LT2E9Z3nc5QaWD7C5HKZXH6Vo3HjZm7Ri+2PAFjnAH1UtV\nt7scpZTKMQ36NEatHsWkzZNw9HPo5t1KqYChQe/y/k/vM23bNBz9HFQpUcXucpRSymM06IF3V77L\njB0zWBG+grtL3G13OUop5VF5OuhFhHdWvsOcX+bgCHdQqXglu0tSSimPy7NBLyKMWDGC73d/z4rw\nFVQsVtHukpRSyivyZNCLCMOjh7No3yKW911OhWIV7C5JKaW8Js8FvYjw5rI3WXpgKdF9oylftLzd\nJSmllFflqaAXEd6IegPHYQfRfaMpW6Ss3SUppZTX5ZmgFxFeW/Iaa46sYVmfZZQpUsbukpRSyify\nRNCLCK8sfoX1x9aztM9SSt9V2u6SlFLKZwI+6J3iZMiiIWw+sZmo3lGUuquU3SUppZRPBXTQO8XJ\niz+8yI5TO1jSewklC5e0uySllPK5gA16pzgZvGAwu87s4sff/0iJwiXsLkkppWwRkEGf4kxh4IKB\nHDh3gMW/X0zxQsXtLkkppWwTcEGf4kyh//z+HL14lEW9FlGsUDG7S1JKKVu5s5UgxpgOxphdxpg9\nxpg3M3i9lzFmq+tjtTGmvudLzVqyM5nw78M5dukYP/T6QUNeKaVwb3PwfMAeoB1wDGsP2R4isivN\nMc2AOBG5YIzpgLWZeLMM2vLaDlPJzmT6zOvDmYQzzO8xnyIFi3jlPEop5Wu+2By8CbBXRA6LSBIw\nE+iU9gARiRGRC66HMUDV7BaUHUkpSfz+u99z/tp5DXmllErHnaCvChxJ8/gotw/ygcDinBR1J5JS\nkug5tyeXEy8zr/s8DXmllErHoxdjjTFtgP5Ay8yOiYiISP08JCSEkJCQbJ8vMSWR7t92J8WZwnfd\nvqNwgcLZbksppfyFw+HA4XB4rD13xuibYY25d3A9HgaIiIxOd9xDwFygg4jsz6Qtj43RX0++Trdv\nu2EwzO46m0L5C3mkXaWU8je+GKOPBeoYY4KMMYWAHkBkuiJqYIV8n8xC3pOuJ1/nudnPkd/k15BX\nSqksZDl0IyIpxpghQBTWfwwTRSTOGDPYelm+Av4OlAU+M8YYIElEmnij4GvJ1+gyqwvFChXjv13+\nS8H8Bb1xGqWUChhZDt149GQ5HLq5mnSVZ2c9S+m7SjO9y3QK5Au4+72UUuo3fDF04xcSkhLoOLMj\n5YqW05BXSqk7kCuC/kriFZ6Z8QyVi1dmauepGvJKKXUH/D7oryRe4ekZT1O9ZHUmd5pM/nz57S5J\nKaVyFb8O+suJl3niP09Qq3QtJnacqCGvlFLZ4LdBf+n6JTpM70C98vX4uuPXGvJKKZVNfhn0F65d\nIGx6GPUr1ueLp78gn/HLMpVSKlfwuwS9EfINKzfks6c+05BXSqkc8qsUPX/tPI9Pe5wmVZvw6ZOf\nYt17pZRSKif8JujPXj1L+6ntaVG9BeM6jNOQV0opD/GLoD+TcIb2U9vTJrgNY8PGasgrpZQH2R70\npxNO025qO0Jrh/LB4x9oyCullIfZGvSnrpyi7ZS2PHXPU/yz3T815JVSygtsC/qTl0/SZkobOtfr\nzHtt39OQV0opL7El6E9cPkGbKW3oen9X3m3zroa8Ukp5kc+D/tilY4RMDqFX/V68HfK2r0+vlFJ5\njltBb4zpYIzZZYzZY4x5M4PX6xpj1hpjrhljht6urZDJIYQ3COdvrf6W3ZqVUkrdgSyD3hiTD/gU\nCAMeAHoaY+qlO+wM8EfgX1m1N/Dhgbz12FvZKNU7PLkBryf5Y11ak3u0Jvf5Y13+WFNOudOjbwLs\nFZHDIpIEzAQ6pT1ARE6LyEYgOavG/tLiL9kq1Fv89Zvqj3VpTe7Rmtznj3X5Y0055U7QVwWOpHl8\n1PWcUkqpXMD2G6aUUkp5V5abgxtjmgERItLB9XgYICIyOoNj3wYuicjYTNry3U7kSikVQHKyObg7\nm6/GAnWMMUHAcaAH0PM2x2daTE4KVUoplT1Z9ujBml4JjMMa6pkoIqOMMYOxevZfGWMqARuAEoAT\nuAzcLyKXvVe6Ukopd7gV9EoppXIvr12MNcZMNMacNMZsS/NcGWNMlDFmtzFmiTGmlLfOn0lN1Ywx\ny40xO40x240xr9hdlzGmsDFmnTFms6uukXbXlKa2fMaYTcaYSH+oyRhzyBiz1fVerfeHmlw1lDLG\nzDHGxLm+h01t/pm61/UebXL9ecEY84rd75Ux5i3X+7PNGPMfY0whP6jpT64ssDUP7jQvXe/lXtfP\nXGhW7Xtz1s03WDdZpTUMWCYidYHlgK/vnEoGhorIA0Bz4GXXzV+21SUi14E2ItIQeAhoa4xpYWdN\nafwJ+CXNY7trcgIhItJQRJr4SU1gDWsuEpH7gAbALjvrEpE9rvfoYeAR4Aowz86aXNf4BgENReQh\nrOuDPW2u6QFgANAI+B3wtDGmtk01uZ2Xxpj7gW7AfcATwGcmqwXDRMRrH0AQsC3N411AJdfnlYFd\n3jy/G/V9D7T3l7qAosB64H67awKqAUuBECDSH75/wEGgXLrn7K6pJLA/g+f95WcqFFhld01AGdf5\ny2CFfKTd//aA54Gv0zz+G/BnIM6OmtzNS6z/AN5Mc9xioOnt2vb1PPqKInISQEROABV9fP5Uxphg\nrP/FY7DeTNvqcg2RbAZOAA4R+cXumoCPsH7o017EsbsmAZYaY2KNMQP9pKaawGljzDeuoZKvjDFF\n/aCuG7oD/3V9bltNInIO+BD4HxAPXBCRZXbWBOwAHnMNkRQFngSq21xTWpnlZfqbWOPJ4iZWu2+Y\nsuVKsDGmOPAt8CexZgalr8OndYmIU6yhm2pYP3ghdtZkjHkKOCkiW7jNdFl8//1rIdZwxJNYw26P\nZVCDr2sqADwMjHfVdgWrx2V3XRhjCgIdgTmZ1ODLn6lawGtYvdYqQDFjzO/trElEdgGjsX5zXQRs\nBlIyOtRXNWUh23X4OuhPuqZiYoypDPzq4/NjjCmAFfLTRGS+v9QFICIXsX7gGtlcUwugozHmADAD\n67rBNOCEne+TiBx3/XkKa9itCfZ/744CR0Rkg+vxXKzgt7susMZvN4rIaddjO2tqBKwRkbMikoJ1\nzeBRm2tCRL4RkUYiEgKcB3bbXVMamdURj/Wbxw3VXM9lyttBb7i1RxgJ9HN9Hg7MT/8FPjAJ+EVE\nxqV5zra6jDHlb1xNN8YUAR7H6lnYVpOIDBeRGiJSC+sGueUi0gdYYFdNxpiirt/EMMYUwxp73o7N\nP1OuX62PGGPudT3VDthpd10uPbH+o77Bzpp2A82MMXe5Lhy2w7rQb+v7ZIyp4PqzBvAs1jCXXTW5\nm5eRQA/XrKWaQB2sa3uZ8+KFhf8Cx4DrWONy/bEuxCzD+qZHAaV9cZEjTU0tsH4124IVppuADkBZ\nu+oC6rvq2AxsBd5wPW9bTenqa83Ni7F2vk8103zftgPD7K4pTW0NsO4g3wJ8B5Syuy6sC/ungBJp\nnrO7pj9j/Se4DZgCFPSDmn7CGqvfjDWjy5b36U7zEmsGzj6sC8ehWbWvN0wppVSAs/tirFJKKS/T\noFdKqQCnQa+UUgFOg14ppQKcBr1SSgU4DXqllApwGvRKKRXgNOiVUirA/T+al80hdsb74wAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1092ff2b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(s_values, s_quantiles, label=\"Source\")\n",
    "plt.plot(r_values, r_quantiles, label=\"Reference\")\n",
    "plt.legend(loc=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To put it simplistically, we need to adjust `source` so that the blue line matches the green line as closely as possible.\n",
    "\n",
    "To do this, we can linearly interpolate on the arrays. This is the **core concept** of the algorithm. The value of the source CDF at each unique input is mapped to the reference CDF (It need not be an exact equality as we use linear interpolation). The reference cdf value is then used to look up the corresponding reference value (interpolated between values if the match isn't exact).\n",
    "\n",
    "The interpolated values have the same shape as `s_values` - the result is the new set of source values to replace the original. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.33333333  0.55555556  0.77777778  0.88888889  1.        ]\n",
      "[ 0.11111111  0.33333333  0.44444444  0.77777778  0.88888889  1.        ]\n",
      "[10 25 35 63 90 99]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 25.        ,  44.33333333,  63.        ,  90.        ,  99.        ])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(s_quantiles)\n",
    "print(r_quantiles)\n",
    "print(r_values)\n",
    "\n",
    "interp_r_values = np.interp(s_quantiles, r_quantiles, r_values)\n",
    "interp_r_values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a way you could think of it as a lookup table mapping from source to target values, independent of it's position in the original array. e.g. `10` becomes `25.0`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{10: 25.0, 20: 44.333333333333336, 40: 63.0, 60: 90.0, 90: 99.0}"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dict(zip(s_values, interp_r_values))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can take our `s_idx` to recreate an array of the same size as the source but with new values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 25.        ,  44.33333333,  63.        ],\n",
       "       [ 44.33333333,  25.        ,  90.        ],\n",
       "       [ 63.        ,  99.        ,  25.        ]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "target = interp_r_values[s_idx].reshape(orig_shape)\n",
    "target"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`target` has the same shape and *relative* spatial pattern as `source` but \n",
    "with values adjusted to match the histogram of `reference`. Let's look at them side-by-side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Source\n",
      " [[10 20 40]\n",
      " [20 10 60]\n",
      " [40 90 10]]\n",
      "Target\n",
      " [[25 44 63]\n",
      " [44 25 90]\n",
      " [63 99 25]]\n"
     ]
    }
   ],
   "source": [
    "print(\"Source\\n\", source2D.astype('uint8'))\n",
    "print(\"Target\\n\", target.astype('uint8'))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also visually assess the quality of the histogram match"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1093f0198>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FOXXxvHvk9ASEgi9BEhoggUUkCaKAWkigiIdAVEU\nG4rY0PenIlYsYMOC0kQFQUQgtNACitKkY+hFCL2TQNruef+YEEIMJITdnc3mfK5rL7ZMZu5swsns\nM08xIoJSSinf5Wd3AKWUUu6lhV4ppXycFnqllPJxWuiVUsrHaaFXSikfp4VeKaV8XJaF3hgz2hhz\n2Biz4QrbfGaM2W6MWWeMucW1EZVSSl2L7JzRjwVaX+5FY8zdQFURqQ70B752UTallFIukGWhF5E/\ngJNX2KQD8H3qtiuAosaYMq6Jp5RS6lq5oo0+FNiX7nFs6nNKKaW8gF6MVUopH5fPBfuIBSqme1wh\n9bn/MMboxDpKKZUDImJy+rXZPaM3qbfMzAB6AxhjGgGnROTw5XYkIl51e+ONN2zPkFtyaSbNlNty\nOZ3Crl3C2LHCQw8J4eFCqVJCp07C558LGzYIDof979XR+KM0H92UNTeVJOH5Z5GUlEtev1ZZntEb\nY34CIoASxph/gTeAAlbNllEiMtsY09YYswOIB/pecyqllMoBEdi5E5Ysgeho69/kZLjzTuv24otw\n/fVgcnxu7Hqbjmyiw6QOdLq+E7VXz8O/YCGXHyPLQi8iPbKxzdOuiaOUUtknAtu2XVrYjblY2P/3\nP7juOu8q7OnN2DqDR2Y8wojWI3iw9oNuO44r2uhztYiICLsjZMobc2mm7NFM2Xe1uUQgJuZiUV+6\nFAoUsIr6XXfBW29BlSrXVtg98V7J7Nl86ljGR9vHE9k9koYVGrr1eMYV7T/ZPpgx4snjKaVyN6cT\nNm2yivqFwh4UZBX2iAjr3/Bwu1NeBRGS33kL/9ffYF21IMr8tZHQEuFZfpkxBrmGi7F5/oxeKeU9\nHA7YsOFiYf/9dyhWzCro990HI0ZAxYpZ78crnT/Pud7dCfxlOk4Dtfq8RP7iYR45tJ7RK6Vsk5IC\n69ZdLOx//AGlS188W2/aFEJ9YfhlbCzxbVtQeMMWEgMKUGDSFEz79tn+8ms9o9dCr5TymORkWLPm\nYmFftgwqVLh48bRpUyhb1u6Urrfu6U7cMnIq8RXLUnjOArjxxqv6ei30SimvlZQEq1df7BXz119Q\nufKlhb1UKbtTuo/D6eDVha8yddNk/oxtQ+lX3oYSJa56P1rolVJeIzERVq68WNhXrIDq1S8W9jvu\nyFGdy5XOJJ6hx9QexCfHM6XzFEoGlszxvrTQK6Vsc/68VcwvdHdcvRpq1rzYK+b22yEkxO6Unrfj\nxA7aT2zPnWF38tndn5HfP/817U8LvVLKY86dgz//vNjGvmYN3HTTxYunTZpAkSJ2p7TRjh2ceKgr\nEc328kT7t3ii/hMu2a0WeqWU28TFWRdMLxT29evh5psvFvbbbrP6tSuQBQtIfKADhc6c49D9LSn7\na5TL9q2FXinlMmfOWF0cLxT2TZugXr2LbeyNG0NgoN0pvYwIKZ9/inluEP5OIb51MwpP/s2lH220\n0CulcuzUKWtQ0oXCvmUL1K9/sbA3bAgBAXan9GIinH/sYQK+GwdA4gsDKfj+R+Dv79LD6MhYpVS2\nHT9+sbBHR8OOHdCokVXUhw+HBg2gYEG7U+YeG49sYn7sbwzI74//6LEU7NXL7kiZ0jN6pXzY0aPW\n/DAXCvuePVa7+oVeMfXqWZOCqas3fct0+s3sxyetRtCzUH2oUcNtx9KmG6VUmsOHL52yNzbW6glz\n4eJp3bqQTz/HXxMR4d3f3+Wr1V8xres06ofWd/sxtdArlYfFxl5sX1+yBI4csfquXyjst9zi8ubi\nvEuE8zu38vD6N9l5Yie/dfuN8sHlPXJoj7TRG2PaAJ9gLT04WkSGZXg9BBgDVAXOAw+LyD85DaWU\nyty+fRfP1pcsgZMnrWkE7rwTnngCatXSwu4W585xrlc3EqPmUPrDtox5ZAkB+XPPVeosz+iNMX7A\nNuAu4ACwCugmIlvSbfMBcFZE3jLG1ABGikiLTPalZ/RKXYU9ey4t7HFxF3vE3HmnNTeWX3ZXflY5\ns38/cW1bErQxdebJGbMwLf5T3tzKE2f0DYDtIrI39YCTgA7AlnTb3AC8ByAiW40x4caYUiJyNKfB\nlMpr0q93euGWlOTd6536vOXLOd+uDUHHT+d45klvkJ1CHwrsS/d4P1bxT2890BFYZoxpAFQCKgBa\n6JW6jPTrnV64waXrnVavroXdLo6DB3BENCUgMZm4Jg0Imj47187I5qrr7+8Dnxpj1gAbgbWAI7MN\nhwwZknY/IiLCa9eyVMrVLqx3mr6wp1/vdOjQa1/vVLnG6YTT9Fj8KPfcX4mHi9xJ0BdfQ/5rm5js\nakRHRxMdHe2y/WWnjb4RMERE2qQ+HgxIxguyGb5mN1BLROIyPK9t9CrPcDph8+ZLC3tQ0MUeMblu\nvdM84sLMk80rN2dEq+Hkz2f/QAO3d680xvgDW7Euxh4EVgLdRSQm3TZFgXMikmyMeRRoIiIPZbIv\nLfTKp+3ZA9OnWxdQL6x3mr6w59r1Tn1ZYiL8+Sfnq1QiMm4NT895mqERQ+l/a3+7k6Vx+8VYEXEY\nY54GorjYvTLGGNPfellGAdcD440xTmAz8EhOAymV25w9C1OmwPffW2fw998PXbrAyJFQ3jPdrNXV\nEIGtW5F58zg78xcKLVtBgYRkXm9XgBUdGzK502TuDL/T7pQupQOmlMoBhwMWLYLx4yEy0jpr79MH\n7rlHpxTwZofiDnHw1Weo8/mUS54/c10Y+Qe9RED/J21KdmU6qZlSHhQTY525//ADlCljFfcRI3x7\n3dNcKTkZDhwgIbQMf/z7B1E7o4jaGcXe03t5unhtaoYE4WjRnKB2HaFlS4r4+EcvPaNXKgvHj8Ok\nSdbZ+/798OCDVoHPhd2pfdvOncjcuZyN/JWCS5exv1QB6vSHWmVq0apKK1pVbUX90PrkI3WEWS4a\naaZz3SjlBsnJMHu2dfa+cCG0bWsV97vu0knBvM3x2B0UbNSEoP1HLnn+TJUKOFevJKRYOZuSuY42\n3SjlIiKwdq115j5xojXrbJ8+MGYMFC1qdzqFwwF+fiQ5k/lr319Wc8yuKLYd38aWpBTyFQkkudmd\nBLXriGnViiKVKtmd2GvoGb3K8w4etNrcv/8e4uOhd2/o1QuqVrU7mWLvXmTePOJm/Ur+xUsZ/Ep9\nxpp11ChRg1ZVreaYxhUakz/2IISG+uyMbtp0o1QOnD9v9XcfPx6WL4eOHa2z99tvz1VNtz7r3PBh\nOD7/jOA9By55/u9X+xL+fx9SIjB3TkWQU9p0o1Q2icCyZdaZ+y+/WGuj9ukDU6fqgtd2S3GmsGL/\nirTmmLZz1vDaniQSCxci8c7bCU5tjqmnH7NyRM/olc/bs8cq7t9/b/Vx79PH6jkTGmp3sjwqNhbm\nzycu8le2FHPwbouCLN6zmPCQcFpVaUXraq1pIhUpePCItTq5Xv3WphulMpNxtGq3blbb+6236qRh\ntti/n8QP3ydx9kyK7Pg37ekDoUVYPPtLWlRpQZmgMjYG9G5a6JVKpaNVvYvD6eDvg38TtTOKtatm\nMvWFlQAkBRQg4fZGBLd7ANOqFdSsaXNS76eFXuV5MTFWcf/hByhb1jpz795dR6t61OHDsGABcVGR\nTH6qGfP+XciCXQsoH1w+bbBSs1/XUKDhbdC4sf7lvUpa6FWelNlo1d694aab7E6WhyxZQtLM6Zyf\nPZ2iMbvSnn7jnZZUbtuDllVaElpEL4S4ghZ6lWdcGK06fvylo1VbtNDrdZ7iFCfrD61n3s55tHj8\nA27dfBKA5IL5OXdbfYLbPYBf9+5QLvePRvUmWuiVT7vcaNXOnXW0qtsdOwYLFnCsbFFmhRwhalcU\n83fOp3hAcVpVbUXfTfm58ZCDAne3swYgFCpkd2KfpYVe+SQdrWqDpCRYvpzk2ZGcmz2d4E3b8BP4\nvn4BZr7YPq2tPSwkzO6keY4WeuUzdLSqPUSETUc2sfebYbR748e051Py+xPfsA7BvR7F77HHbEyo\nPDIy1hjTBviEiytMDcvwegngB6Ac4A98LCLjchpK5R0XRquOH2+NUL31Vh2t6lYnT8LmzRypcx0L\ndi1Im6e9UL5CdKzalDuqV6JQm3YUbHsv+Zo2paj+EHxCdtaM9QO2Ya0ZewBYBXQTkS3ptnkDKCQi\nrxhjSmKtMVtGRFIy7EvP6BWgo1U9JiUFVq4kZe5s4mdNI3hdDIn+EPZ6MLdVa0brqq1pVbUVVYtr\nm5g388QZfQNgu4jsTT3gJKADsCXdNoeAWqn3g4HjGYu8UmfOWHPMXBit2rWrdYFVR6u6noiw9WgM\nobWaEHzkFPmAooDD34+kurWI7TmV/JW1uOcV2Sn0ocC+dI/3YxX/9L4FFhpjDgBBQFfXxFO5XWaj\nVZ991uoaWbCg3el8xNmz4O/PCZNwSXMMwIzwIlQNDCB/m3so1LY9/hERFA0Otjmw8jRX9T5+BVgv\nIs2MMVWB+caY2iISl3HDIUOGpN2PiIggIiLCRRGUN8lstKqureoiDgesWYNj7hziIqcS9Pcm3uoT\nzvBqR2ka1pRWVVvxwm0vUKNEDczjCRAQYHdidZWio6OJjo522f6y00bfCBgiIm1SHw8GJP0FWWPM\nbOAdEVmW+ngh8LKIrM6wL22j92E6WtX9Dn/2HkVfe4dCZ+LTnnP6Gf4d0JtyH31DwXz6MckXeaKN\nfhVQzRgTBhwEugHdM2wTA7QAlhljygDXAbtQPi8pCebMuXS06tChOlr1momAMZxOOM2i3YvS5mm/\nc+1JxpyJJ658Kfzb3E3APR3wa96c8JAQuxMrL5atfvSp3Ss/5WL3yveNMf2xzuxHpfa0GQtUAgzw\nnohMzGQ/ekbvA0RgzRrroqqOVnURpxM2bEhrjjmafJo+/Uuz4fAGbqt4W1rvmBsDwzGHDlkjx/QK\ndp6hA6aUx1wYrTp+/MXRqr1762jVa3LyJHH9H8Zv4UICT5xNezq5QD6Wrv2N26o3JyC/trHndVro\nldudOAFPPgnz5uloVVdauGshz8x6it8Hb6P4OSG+THFM69YE3nMf3HUXlMhb66Kqy9M1Y5Vb/fmn\nNbf7Aw9YF1gLF7Y7Ue4nInz818d8/NfHjGk/hpAqDqhSlcI1a2pzjHILLfQqU04nDBsGn3wC330H\n995rdyLfEJ8UzyMzHmHHiR2s6LeCSkUrQXW7Uylfp4Ve/cfhw1bb+7lzsHo1VKxodyIfcOoUp155\njlY1VnJj5Qb83vd3bXtXHqOtrOoSixZB3bpQvz4sXqxF3iU2bSL+5hsI+XocY5cWZ0z7MVrklUfp\nxVgFWHNfDR1qNdN8/73VD15dO+fPk0h5qDcFEpKJu6E6QZHzoHJlu2OpXEYvxqprFhsLPXpA/vxW\n//iyZe1O5AOcThJffI6Cwz+jAHCua0eCxkzQuZeVLbTpJo+bNQvq1YNWrazuk1rkXSPm2Bbm/DEO\nh58hZcTHBE78RYu8so2e0edRSUnwf/8HP/8MU6bAHXfYnch3TIuZRv/I/nz4yTD8zS3QqJHdkVQe\np4U+D9q9G7p1g9KlraaakiXtTuQbHE4Hry9+nQkbJjCrxyzqh9a3O5JSgDbd5DlTp0LDhlahnzFD\ni7xLJCdzasdm2k1sx5/7/2T1Y6u1yCuvooU+j0hIgKeegpdestrln3tOB2G6xKFDxDVtxLHb63BL\nYFXm95pP6cKl7U6l1CW00OcB27ZZzcRHj1pNNfX1ZNM1li/n3M03ELR8DeUliPeue4J8ftoaqryP\nFnof98MP0KQJPPGEdeFVpxF2DcfXX5HS9HYCj5wkrmFdAtf/AzfeaHcspTKlA6Z8VHw8DBhgTUr2\n889w8812J/Idp+bNIKRNBwASnniUQp+OtAYhKOUm1zpgSs/ofdDGjVbzjNNpzVWjRd51/j7wNzfH\nDODP+2/FOXYshb4cpUVeeb1sFXpjTBtjzBZjzDZjzMuZvP6CMWatMWaNMWajMSbFGKNrm3mYCHz7\nLTRvDi+/DOPGQVCQ3al8x7h142jzYxtGtB7Bbb+uwu+hh+yOpFS2ZGdxcD9gG3AXcABrDdluIrLl\nMtu3AwaKyH9mS9GmG/c5cwYeewxiYqymmpo17U7kO5IcSTw39zkW7F7AtK7TuKHUDXZHUnmMJ5pu\nGgDbRWSviCQDk4AOV9i+O/Cf9WKV+/z9tzXjZLFisHy5FnmXOX+e+Id78cg7Ddl/dj8r+63UIq9y\npez0BQsF9qV7vB+r+P+HMSYAaAM8de3RVFZE4LPP4J13YORIa3Fu5SJ79hB3b2uCNm3jw6plKb11\nFX7+2nVS5U6u/s29F/hDRE5dboMhQ4ak3Y+IiCAiIsLFEfKGEyegb184cMA6i69Sxe5EvkPmzyeh\n8/0EnY4nvlI5yk6bB1rklQdFR0cTHR3tsv1lp42+ETBERNqkPh4MiIgMy2TbX4HJIjLpMvvSNnoX\nuLCOa6dO8N57UKCA3Yl8R/Lwj/B/8SX8nEJ8i6YUnvyb1SamlI080Ua/CqhmjAkzxhQAugEzMglS\nFLgTmJ7TMOrKnE6rsN9/P3zxBXz8sRZ5V9p3eh/vbPoSP6eQ9MpLFJ67SIu88glZfh4VEYcx5mkg\nCusPw2gRiTHG9LdellGpm94HzBOR8+6Lm3fpOq7uFb0nmu5TuzPokUHIs60ooIMPlA/RkbG5wKJF\n0KuX1SY/ZAjk0+ZilxERPln+CcOWDWPC/RNoWbWl3ZGU+g9dStCH6TqubuR0cn7tKvrt/YyYozEs\n77ec8JBwu1Mp5RY6BYKXio2Fu+6yLryuWaNF3qVOnya+XWtM48ZU2X2aZQ8v0yKvfJoWei+k67i6\nUUwM8XVuovCcBVCoEENvHkhA/gC7U3lMeHg4xhi9eektPDzcLT93baP3IunXcf3pJ7j9drsT+RaZ\nOpWkXj0oeD6JuJpVCIqMgqpV7Y7lUaltvXbHUJdxuZ+PttH7iPTruK5dCyVK2J3It8Tt341/z64E\nJDo490B7gsb/BIUL2x1LKY/QphsvkHEdVy3yrrXt+DYazmjHhGfuJOWD9wmc8psWeZWnaNONjRIS\n4PnnYe5cmDRJl/hzhxlbZ9BvRj/eaf4Oj9Z71O44ttOmG++mTTc+Zts26NIFrrvO6lWjS/y5llOc\nvBn9JmPWjWFG9xk0qtDI7khK2Uabbmyg67i6UXIyCQOe5IvHbmHxnsWsfnS1Fvlc5I8//qBJkyaE\nhIRQsmRJ7rjjDv7++2+7Y+V62nTjQenXcZ08GWrXtjuRjzlyhPj776Hwn6tJCMiP/55/yV9a+6am\n581NN2fPnqVSpUp88803dO7cmaSkJH7//XfKli3LTTfd5LLjOBwO/P39XbY/V3JX042e0XtIxnVc\ntci72KpVnKt9A4X/XM25kkUpFLVIi3wus23bNowxdOnSBWMMBQsWpEWLFtx0002ICG+//Tbh4eGU\nLVuWhx56iLNnzwKwZMkSKmaY/Kly5cosWrQIgDfffJPOnTvTq1cvQkJCGD9+PE6nk3fffZdq1apR\ntGhR6tevT2xsLABbtmyhVatWlChRguuvv54pU6Z49o1wBxHx2M06XN7idIqMGiVSsqTIuHF2p/FN\nKZMnS1J+fxGQs/Vqi8TG2h3Ja3nz/8EzZ85IyZIlpU+fPjJnzhw5efJk2mujR4+W6tWry549eyQ+\nPl46duwovXr1EhGR6OhoqVix4iX7Cg8Pl4ULF4qIyJAhQ6RAgQIyY8YMERFJSEiQDz74QGrXri3b\nt28XEZENGzbIiRMnJD4+XipWrCjjx48Xp9Mp69atk1KlSklMTIwn3oLL/nxSn89x7dUzejc6c8aa\nN/6LL+D336FPH7sT+Z7j547z6K4RJPnD+X4PEfTnKihf3u5YuZYxrrnlRHBwMH/88Qd+fn489thj\nlCpVivvuu48jR47w008/MWjQIMLCwggMDOS9997j559/xul0ZmvfjRs35t577wWgYMGCjB49mnfe\neYdq1aoBUKtWLYoVK0ZkZCSVK1emd+/eGGO4+eab6dixY64/q9dC7ya6jqv7rT24llu/vZWSdW+n\nYMw2Ar4dqxP0XyMR19xyqkaNGowZM4Z///2XzZs3c+DAAQYOHMjBgwcJCwtL2y4sLIzk5GQOHz6c\nrf1mbNrZt28fVTJZlm3v3r0sX76c4sWLU7x4cYoVK8ZPP/3EoUOHcv5NeQHtXuliuo6rZ0xYP4FB\nUYMY2XYkXW7sYncc5QbXXXcdffr04ZtvvqF8+fLs3bs37bW9e/eSP39+ypQpQ2xsLOfOnUt7zeFw\ncPTo0Uv2ZTJ8zKhUqRI7d+7khhsuXey9YsWKREREMG/ePDd8R/bJ1hm9MaaNMWaLMWabMebly2wT\nYYxZa4zZZIxZ7NqYucOJE3DffVb3yeXLtci7nAgp8+by7OxnGLp0KIv7LNYi70O2bt3K8OHD0y6K\n7tu3j4kTJ9K4cWO6d+/OiBEj2LNnD3Fxcfzf//0f3bp1w8/Pj+uuu46EhATmzJlDSkoKb7/9NklJ\nSVc81iOPPMJrr73Gjh07ANi4cSMnT56kXbt2bNu2jR9++IGUlBSSk5NZvXo1W7Zscfv3705ZFnpj\njB/wBdAauBHoboypmWGbosBIoJ2I3ATkuRL3559Qpw5UqwbLluli3S6XkMD53j3I1+ZuakxeyKpH\nV3FTadd1uVP2Cw4OZsWKFTRs2JDg4GBuu+02ateuzccff0zfvn3p1asXTZs2pWrVqgQGBvLZZ58B\nUKRIEb788kseeeQRKlSoQHBwMBUqVLjisQYNGkSXLl1o1aoVRYsWpV+/fpw/f56goCCioqKYNGkS\n5cuXp3z58gwePDjLPxzeLruLg78hInenPv7P4uDGmCeAciLyehb7kqyOl9s4nTBsGHz6qbVASLt2\ndifyQf/+S9y9bQjaEENSwXzkGzMOvx497U6VK3lzP3pl7xQIocC+dI/3Aw0ybHMdkD+1ySYI+ExE\nJuQ0VG6Rfh3XVat0HVe3WLyY8x3bE3QqjvgKZSgcOQ90PVelroqret3kA+oCdwNtgNeMMdVctG+v\ntGiR1aumfn1YvFiLvDskJp3n336dCTgVR3xEEwqv/0eLvFI5kJ0z+ligUrrHFVKfS28/cExEEoAE\nY8xS4GZgR8adDRkyJO1+REQEERERV5fYZrqOq2fEnoml05RO1Hm2HiOO1KXwm2+Dlw5bV8rVoqOj\niY6Odtn+stNG7w9sBe4CDgIrge4iEpNum5rA51hn8wWBFUBXEfknw75ydRt9bCz06GF11Z4wQZf4\nc5ff9/5O11+6MqDBAAbfPvg/XeNUzmkbvXezba4bEXEATwNRwGZgkojEGGP6G2MeS91mCzAP2AAs\nB0ZlLPK53YV1XFu31nVc3UWcTj5f8TmdpnRiTIcxvHLHK1rklXIBnb0yC7qOqweIkPzOW6yaP44n\nugUxrdtvVCmm/VPdQc/ovZsuPGIDXcfVA86eJb5nFwrPnMttwPK3FhKgRV4pl9K5bi5D13H1gK1b\nia97E4VnziUhqBAyYwYBTZvbnUopn6OFPoOEBHjqKXjpJatd/rnncj4bn7o8WbWKxHq3UHjHv8RX\nD6PQmg2Y1NkFlbpaX331FWXLlqVIkSKcPHnS7jheRwt9Otu2QaNGcPSotY6rLtbtHnFJcfTa9j67\nSvhxrkNbCq/ZBNWr2x1L2Sw8PJzAwECKFClCuXLl6NWrV9riIleSkpLC888/z4IFCzhz5gzFihXz\nQNrcRQt9Kl3H1TN2nNhB49GNKRBUlMqrdxA4LRKCguyOpbyAMYZZs2Zx5swZ1q9fz8aNG3n77bez\n/LpDhw6RmJjI9ddfn6PjZndO+9wszxf6+Hh4+GF4+21YuBD699emGneZvX02TcY04clbn2R0+9EU\nKlVO32x1iQs9TkqXLk3r1q3ZvHkzAElJSbzwwguEhYVRrlw5nnzySRITE9m+fTs1Uxd7KFasGC1S\nRzBeaTnAvn378uSTT3LPPfcQHBxMdHT0ZfcPF5cqHD58OGXKlCE0NJRx48al7S8hIYHnn3+e8PBw\nihUrRtOmTdO+dvny5TRp0oRixYpRp04dlixZ4vb3MFPXsjzV1d7wsmXMjh4Vuf56kT59RM6etTuN\n73JMmSzDIl+V0I9DZdm/y+yOk6d52//B9NIv/7dv3z6pVauWDB06VEREBg4cKB06dJBTp05JXFyc\ntG/fXl599VUREdmzZ4/4+fmJ0+kUEcl0OcCSJUumLQf40EMPSUhIiPz1118iYi0teKX9R0dHS758\n+WTIkCGSkpIis2fPlsDAQDl16pSIiDz55JPSrFkzOXjwoDidTvnrr78kKSlJYmNjpUSJEjJ37lwR\nEVmwYIGUKFFCjh07dtn34HI/H65xKcE83Y++SxeoUAGGD7c7iY9KTCTxmacoOGo00XWLU2PxRsoV\n0WX+7JRVP3rzpms+YckbV///vHLlyhw/fhyAuLg4OnTowNSpU/Hz8yMoKIiNGzdSuXJlAP766y96\n9uzJrl272LNnD1WrViU5ORk/Pz8mT57MyJEjLzl7fvzxxwkNDeW1116jb9++iMglZ+VX2v+SJUto\n27YtZ8+exc/PagQpU6YMM2fOpH79+hQuXJiVK1dy002XTpv9wQcfsHnzZsaPH5/2XJs2bejZsye9\nevXK9D3QfvQu9vPPsHGjNV+NcoM9ezh/XzsC1m8mOZ8ft/cbSr7gcnanUlnISYF2penTp9OsWTOW\nLl1K+/bt+fvvvwkPD+fcuXPUq1cvbTun05lWEDOOnk6/HCBYrRYOh4PevXunbZN+acGjR49ecf8A\nJUqUSCvyAIGBgcTFxXHs2DESExMvuyzh5MmTmTlzZlqOlJQUmjf3fBfiPFnoDx2CZ56BmTOhUCG7\n0/igyEiSenYn4EwcZ8uXJHjGXGv+CKWycKG4Nm3alKeffpqXX36ZhQsXEhgYyObNmylXLuuThews\nB5j+j0NWN+OBAAAagUlEQVTJkiWvav/plSxZkkKFCrFz505q1ar1nxy9e/fmm2++uap9ukOeuxgr\nYl1w7dcPGmScVV9dM4fTwYpRb1DgTBynWtxO8MatWuRVjgwcOJCVK1eycuVKHn30UQYOHJi2Fmxs\nbCxRUVFp26Y/+77ccoBbt27N9DjGmCz3fznGGPr27cugQYM4ePAgTqeT5cuXk5yczIMPPsjMmTOJ\niorC6XSSkJDAkiVLOHDgwLW8LTmS5wr9Dz9YUxu8fsW1sFROnDh/gnt+uofX2gdxZuRwQuYtgdSP\nz0plJWMTTMmSJenTpw/Dhg1j2LBhVKtWjUaNGhESEkKrVq3Ytm1bpl97ueUAL/SEyUxW+79S1o8+\n+ohatWpRv359SpQoweDBg3E6nVSoUIHp06fz7rvvUqpUKcLCwvjoo49s6c6Zpy7GxsZa67rOm2f9\nq1xn/aH1dJzckQ41OvBByw/I55cnWwW9nk5q5t30Yuw1ErGaa556Sou8y4jA4cNMPLqYZ+Y+w6dt\nPqVHrR52p1JKZZBnCv3o0XDkCLz6qt1JfMSpUzj79OHEqiUMeyaEBQ8t4OayusyfUt4oTzTd7N0L\nt95qre2aoauryom1a3E8cD/+u/cSF5gP56xIikS0tjuVygZtuvFutq0wlXqQNsaYLcaYbcaYlzN5\n/U5jzCljzJrU2/9yGsjVnE5rioPnn9cif81E4LvvcDZuhP/uvcRWL0vA+n+0yCvl5bJsujHG+AFf\nYK0ZewBYZYyZLtbygektFZH2bsh4Tb7+2prP5oUX7E7iA1asgEcfxQ/Y3bkllb+foQMRlMoFstNG\n3wDYLiJ7AYwxk4AOQMZC73WzU+3caXWj/OMPyJdnrka4R5IjieeOT6DOXcW4+/6XqPzUYLsjKaWy\nKTvlLxTYl+7xfqzin1FjY8w6IBZ4UWxeHNzphL59rYuvqZPbqRw6ePYgnad0pkRgCd6N3E3RQjqH\ns1K5iasGTP0NVBKRW7CaeX5z0X5z7NNPrSblZ5+1O0kulXpB6M99f1L/2/q0rtqaaV2naZFXKhfK\nzhl9LFAp3eMKqc+lEZG4dPfnGGO+NMYUF5ETGXc2ZMiQtPsRERFERERcZeSsbd0K77xjNSn7+7t8\n975v/36kd29+eaghTx0ezdgOY7nnunvsTqVUnhEdHU10dLTL9pdl90pjjD+wFeti7EFgJdBdRGLS\nbVNGRA6n3m8ATBaR8Ez25fbulQ4H3H47PPigNThKXaWoKKRnT8yxY6yqFkjI8nVUL6HL/PkKb+1e\nGRwcnDatQHx8PAULFsTf3x9jDN988w3du3f3WJbExEQCAgLYv38/5ct7dlpt20bGiojDGPM0EIXV\n1DNaRGKMMf2tl2UU0MkY8wSQDJwHuuY00LX66CMICLCWBFRXweGAt99G3nwTI8K6m8tw/ezlBJUI\ntzuZygPSrw1bpUoVRo8eTbNmzXK0L4fDgf81fJQXkf/Mu5PrXcuqJVd7w82r22zcKFKypMju3W49\njG/q1EkExGGQ3x9pKc7kZLsTKTdw9/9BV0i/0tQFy5Ytk4YNG0pISIiEhobKc889Jw6HQ0SsFaKM\nMfLVV19J1apV5frrrxcRkcjISKlevboUK1ZMBg4cKI0aNZIff/wxbZ9ff/211KhRQ0qUKCHt2rWT\nAwcOiIhIgwYNxM/PTwoXLizBwcEyffp0D33n7lthymdmr0xOhj594N13ITzc7jS5i4gQeUsgxwsb\n1o19j9u/i8Jof1TlRQoUKMDIkSM5efIkv//+O5GRkXz33XeXbDNr1izWrFnD2rVrOXToEN26dePT\nTz/l6NGjlC9fnjVr1qRt+/PPP/PZZ58xa9YsDh8+TJ06dejZsycAS5cuRUTYvn07Z86coX17rxse\ndPWu5a/E1d5w49nE0KEirVuLpC4bqbIpPileekztIXW+riN79qy3O45ysyz/D1r9rf57u5rtr1Fm\nZ/QZvf/++9KjRw8RuXhGv3z58rTXR40aJc2bN0977HQ6pXTp0mln9M2aNZOffvop7fWkpCTJnz+/\nHDlyJG1/sbGx1/y9XK3L/XzQM3pYtw4+/xy++w58rWnNnXad3EXj0Y3xN/4se3gZYWG17Y6kVKZi\nYmJo27YtZcuWpWjRorz11lscO3bskm0qVKiQdv/AgQOXLBdojCE0NDTt8d69e3n88ccpXrw4xYsX\np3Tp0hQoUID9+/e7/5uxQa4v9ElJVpPNhx9aC32rLGzcCBMnMm/HPBqPbky/Ov0Yf994AvIH2J1M\neYPLndNfzfZu8Oijj1KvXj12797N6dOnee211/7TOyX9BdRy5cqxb9/FcZ4iQmzsxV7hFStWZNy4\ncZw4cYITJ05w8uRJ4uLiqFOnju9diMUHCv3QoRAWBunW/VWXM2EC0rAhKX16MXzkg0zpPIUBDQf4\n5C+28i1xcXEULVqUgIAANm/ezLfffnvF7du3b8/KlSuZO3cuDoeDjz/+mFOnTqW9/vjjj/PWW2+l\nrSJ18uRJfv31V8C6HhASEsKuXbvc9w15WK4u9KtWwbffwqhR2mRzRQkJ1kK5vXtjzp9n3q3FGPvy\nXzQNa2p3MqX+I7MTjxEjRvDtt99SpEgRBgwYQLdu3a74NWXLlmXixIkMGDCAUqVKceDAAWrVqkXB\nggUB6NatGwMGDKBjx46EhIRQt25dFixYkPb1Q4cOpVOnThQvXpzIyEg3fJeelWvno09IgLp1rUnL\nMvzMVXq7d0OnTrBmDYn5DD8/fjtdh0dRML/OOpkXeeuAKXdzOByULVuWyMhIGjZsaHecy7J1Pnpv\n9PrrcMMN0NW2oVm5REoKyVtj2FPcj1nj/o/eny/VIq/yhLlz53LmzBkSEhIYMmQIhQsXpl69enbH\nskWu7Cz9558wYQJs2KBNNlfiFCdv7v+BjX2DeOXx6XS8saXdkZTymKVLl9KzZ08cDgc33XQT06ZN\nI18eHR+S65puzp2DW26B99+Hjh1dFMwHHY47TO/fenM++TyTO0+mbFBZuyMpL5BXm25yC226SfXK\nK1C/vhb5TG3aBCIs3LWQuqPqUr98fRb1WaRFXqk8Lld9jlmyBH75xeoKrtJxOuHDD5FXX2XWU614\nLGw939//PS2qtLA7mVLKC+SaQh8XZ60Y9c03ULy43Wm8yMmT1oixmTMxwIl9W1nz7ho9i1dKpck1\nbfRPPGF1qRw71sWhcrO//7a6Tu7Zw8kAQ9SQB+n84jj8TK5rkVMeom303s22+ei9wfz5MGuWNtlc\nQgTn00/jt2cP6yvmJ2nij3Rt0tnuVMrLhYWF6UhoLxYWFuaW/Xp9oT99Gh55xJqwrKguV5pm96k9\nPH/vOXoXq8wdP/5BiWKeXQlH5U579uyxO4KyQbY+4xtj2hhjthhjthljXr7CdvWNMcnGGJf1iRk0\nCO6+G1q1ctUec7+p/0yl4XcNadq8Lx1m7dQir5S6oizP6I0xfsAXWGvGHgBWGWOmi8iWTLZ7H5jn\nqnCzZsGiRdbAqDwvJYUEUnh+3vPM3TmXWT1mUT+0vt2plFK5QHaabhoA20VkL4AxZhLQAdiSYbsB\nwC+AS6rPiRPWPFwTJkBwsCv2mEslJsILL3D23+3c0eYg15WswZrH1lC0kLZjKaWyJzuFPhTYl+7x\nfqzin8YYUx64T0SaGWMueS2nnnnGGhSVw/WBfcPevdClC6xcSUF/+L+u/6NTp6F6MU0pdVVcdTH2\nEyB92/1lK9GQIUPS7kdERBAREfGfbaZNgxUrrJWj8qy5c5GePTEnTnCgWD7ifxxH57t72p1KKeUB\n0dHRREdHu2x/WfajN8Y0AoaISJvUx4Ox1i8clm6bCzP0G6AkEA88JiIzMuwry370R49C7drWCNgm\nTa722/ERkZFI+/YYEdbXDaVq5DKCyrmn25VSyvt5oh/9KqCaMSYMOAh0A7qn30BEqqQLNBaYmbHI\nZ9dTT0HPnnm3yIsIY4v/yy0V/DHt21PnsyngpwOglFI5l2WhFxGHMeZpIAqrO+ZoEYkxxvS3XpZR\nGb8kp2F+/tkaFDV+fE73kLudSTzDYzMfI+ZYDD+vWEPNcrXsjqSU8gFeMwXCoUNw880wcyY0cMnl\n3Nxl9YHVdPulGy2rtGR46+G6WLdSKs21Nt14RaEXgfvvt1aMevddj8Wx3+nTyP/+x9cdQnlj7XBG\nth1J5xt1GgOl1KV8Yq6bH36AXbuspps8Y/16HA90xH/nLsJWl2D57JVUKVYl669TSqmrZPsZfWws\n1KkD8+ZZ/+YJ48bheLw//olJxFYpRanZ0RSocYPdqZRSXipXrzAlAv36WT1t8kSRdzqRRx6Bvn3x\nT0zi3wdaErpprxZ5pZRb2VroR4+Gw4fh1VftTOE5h88dZebe+STm9+PEFx9S6ZcoCNCLrkop97Kt\n6WbvXrj1VmvSslp5oBfhot2L6DWtF4/e0Jv/hXYjX62b7Y6klMolcuXFWKfTmmN+0CDfL/IpzhSG\nLhnKd2u+03VclVK2sKXQf/01nD0LL75ox9E9JDaWw/u20GXLUAr4F2BNf13HVSllD4833ezYITRs\nCH/8ATVreuzQnrVwIYldO3FIzjJ1wisMvPtNXcdVKZVjua7XTd++8MorPlrknU5S3hqKs1VLCh4/\nRdHrb2FQ/We0yCulbOXxCuR0wsCBnj6qBxw/zrnWzcn3+hv4OYVzg58nZMkKKFXK7mRKqTzO4230\n48aBv7+nj+p+f08aQb0FSzhfJJBCE6cQ2Lat3ZGUUgrwgpGxvuC3Lb/RP7I/K871JLzPsxCmc8cr\npVwnV3av9CVT/5nKU7OfYk7POYSXq2t3HKWU+g8t9NdgyuYpDJgzgLkPzuWWsrfYHUcppTKl3UFy\naNKmSTwz9xmiekVpkVdKebVsFXpjTBtjzBZjzDZjzMuZvN7eGLPeGLPWGLPaGNPc9VG9x48bfmTQ\nvEHM7zWf2mVq2x1HKaWuKDuLg/sB24C7gANYa8h2E5Et6bYJFJFzqfdrAdNEpFom+8r1F2MnrJ/A\n4IWDmd9rPjeU0lknlVLu54kBUw2A7SKyV0SSgUlAh/QbXCjyqYKAYzkN5M3GrRvHKwtfYWHvhVrk\nlVK5RnYKfSiwL93j/anPXcIYc58xJgaYDTzjmnjeY8zaMby2+DUW9l5IzZK+OKxXKeWrXNbrRkR+\nA34zxtwOTABqZLbdkCFD0u5HREQQERHhqghuM+rvUby99G0W9V5E9RLV7Y6jlPJx0dHRREdHu2x/\n2WmjbwQMEZE2qY8HAyIiw67wNTuBBiJyPMPzua6N/qtVX/H+svdZ2Hsh1Yr/57KDUkq5nSfa6FcB\n1YwxYcaYAkA3YEaGEFXT3a8LkLHI50YjV45k2LJhLO6zWIu8UirXyrLpRkQcxpingSisPwyjRSTG\nGNPfellGAQ8YY3oDSUA80NWdoT3hsxWf8cnyT4h+KJrwkHC74yilVI7pXDeZGPHXCL5Y9QWLei8i\nLETnrVFK2UvnunGxj/78iK9Xf010n2gqFq1odxyllLpmWujTef+P9xmzdgzRD0VToUgFu+MopZRL\naKFP9c7Sd5iwYQLRD0VTPri83XGUUspltNADQ5cMZeKmiSzus5hyweXsjqOUUi6Vpwu9iPDmkjeZ\n8s8UovtEUyaojN2RlFLK5fJsoRcRXl/8Or9t/Y3FfRZTunBpuyMppZRb5MlCLyK8uvBVZu+YzaLe\niyhVWBfwVkr5rjxX6EWElxe8zPxd81nYeyElA0vaHUkppdwqTxV6EeGFqBeI3hvNwt4LKR5Q3O5I\nSinldnmm0IsIz817jmX7lrGg1wKKBRSzO5JSSnlEnij0IsIzc55h5YGVzO81n5BCIXZHUkopj/H5\nQu8UJ0/Pfpq1h9YS9WAURQsVtTuSUkp5lE8Xeqc4eSLyCTYd3cS8B+dRpGARuyMppZTH+Wyhd4qT\n/jP7s+X4Fub2nEtwwWC7IymllC18stA7nA76zezHrpO7mNNzDkEFguyOpJRStvG5Qu9wOug7vS/7\nz+xndo/ZFC5Q2O5ISillq+wsJYgxpo0xZosxZpsx5uVMXu9hjFmfevvDGFPL9VGzluJMoc9vfThw\n9gCRPSK1yCulFNlbHNwP2AbcBRzAWkO2m4hsSbdNIyBGRE4bY9pgLSbeKJN9uW2FqRRnCr2m9eL4\nueNM7zadgPwBbjmOUkp5micWB28AbBeRvSKSDEwCOqTfQESWi8jp1IfLgdCcBsqJZEcyPX/tyamE\nU1rklVIqg+wU+lBgX7rH+7lyIe8HzLmWUFcj2ZFM96ndiUuKY1rXaVrklVIqA5dejDXGNAP6Ardf\nbpshQ4ak3Y+IiCAiIiLHx0tyJNH1l644nA5+7fIrBfMVzPG+lFLKW0RHRxMdHe2y/WWnjb4RVpt7\nm9THgwERkWEZtqsNTAXaiMjOy+zLZW30iSmJdPmlCwbD5M6TKeBfwCX7VUopb+OJNvpVQDVjTJgx\npgDQDZiRIUQlrCLf63JF3pUSUxJ5YPID+Bt/LfJKKZWFLJtuRMRhjHkaiML6wzBaRGKMMf2tl2UU\n8BpQHPjSGGOAZBFp4I7ACSkJdPy5I4ULFOanjj+R3z+/Ow6jlFI+I8umG5ce7Bqbbs4nn+f+n+8n\npFAIP3T8gXx+PjfeSyml/sMTTTde4VzyOdpPak+JwBJa5JVS6irkikIfnxTPvRPvpWxQWb6/73st\n8kopdRW8vtDHJ8XTbmI7KhapyLgO4/D387c7klJK5SpeXejjkuK4+8e7qRJShdHtR2uRV0qpHPDa\nQn828SxtfmhDzZI1+bb9t1rklVIqh7yy0J9OOE3rH1pTq3Qtvm73NX7GK2MqpVSu4HUV9EKRr1O2\nDl/e86UWeaWUukZeVUVPJZyi5YSWNAhtwBdtv8Aae6WUUupaeE2hP3H+BC2+b0GTik34tM2nWuSV\nUspFvKLQHz93nBbft6BZeDOGtx6uRV4ppVzI9kJ/7Nwx7vr+LlpVbcUHLT/QIq+UUi5ma6E/Gn+U\n5uObc0/1e3jvrve0yCullBvYVugPxx2m2fhm3FfzPt5u/rYWeaWUchNbCv2huEM0G9+Mzjd0Zmiz\noVrklVLKjTxe6A+cPUDEuAh61OrBGxFvePrwSimV52Sr0Btj2hhjthhjthljXs7k9RrGmD+NMQnG\nmEFX2lfEuAj63NyH/zX9X04zK6WUugpZFnpjjB/wBdAauBHoboypmWGz48AA4MOs9tevbj9eueOV\nHER1D1cuwOtK3phLM2WPZso+b8zljZmuVXbO6BsA20Vkr4gkA5OADuk3EJFjIvI3kJLVzl5q8lKO\ngrqLt/5QvTGXZsoezZR93pjLGzNdq+wU+lBgX7rH+1OfU0oplQvYPmBKKaWUe2W5OLgxphEwRETa\npD4eDIiIDMtk2zeAsyIy/DL78txK5Eop5UOuZXHw7Cy+ugqoZowJAw4C3YDuV9j+smGuJahSSqmc\nyfKMHqzulcCnWE09o0XkfWNMf6wz+1HGmDLAaiAYcAJxwA0iEue+6EoppbIjW4VeKaVU7uW2i7HG\nmNHGmMPGmA3pnitmjIkyxmw1xswzxhR11/Evk6mCMWaRMWazMWajMeYZu3MZYwoaY1YYY9am5nrX\n7kzpsvkZY9YYY2Z4QyZjzB5jzPrU92qlN2RKzVDUGDPFGBOT+jNsaPPv1HWp79Ga1H9PG2Oesfu9\nMsa8kvr+bDDG/GiMKeAFmZ5NrQW21oOrrZep7+X21N+5Vlnt3529bsZiDbJKbzCwQERqAIsAT4+c\nSgEGiciNQGPgqdTBX7blEpFEoJmI1AFqA82NMU3szJTOs8A/6R7bnckJRIhIHRFp4CWZwGrWnC0i\n1wM3A1vszCUi21Lfo7pAPSAemGZnptRrfI8CdUSkNtb1we42Z7oReAS4FbgFaGeMqWpTpmzXS2PM\nDUAX4HrgbuBLk9WEYSLithsQBmxI93gLUCb1fllgizuPn418vwEtvCUXEAisBG6wOxNQAZgPRAAz\nvOHnB+wGSmR4zu5MRYCdmTzvLb9TrYDf7c4EFEs9fjGsIj/D7v97QCfg23SP/we8CMTYkSm79RLr\nD8DL6babAzS80r493Y++tIgcBhCRQ0BpDx8/jTEmHOuv+HKsN9O2XKlNJGuBQ0C0iPxjdyZgBNYv\nffqLOHZnEmC+MWaVMaafl2SqDBwzxoxNbSoZZYwJ9IJcF3QFfkq9b1smETkJfAz8C8QCp0VkgZ2Z\ngE3AHalNJIFAW6CizZnSu1y9zDiINZYsBrHaPWDKlivBxpgg4BfgWbF6BmXM4dFcIuIUq+mmAtYv\nXoSdmYwx9wCHRWQdV+gui+d/fk3Eao5oi9XsdkcmGTydKR9QFxiZmi0e64zL7lwYY/ID7YEpl8ng\nyd+pKsBzWGet5YHCxpiedmYSkS3AMKxPrrOBtYAjs009lSkLOc7h6UJ/OLUrJsaYssARDx8fY0w+\nrCI/QUSme0suABE5g/ULd6vNmZoA7Y0xu4CJWNcNJgCH7HyfRORg6r9HsZrdGmD/z24/sE9EVqc+\nnopV+O3OBVb77d8iciz1sZ2ZbgWWicgJEXFgXTO4zeZMiMhYEblVRCKAU8BWuzOlc7kcsVifPC6o\nkPrcZbm70BsuPSOcATyUer8PMD3jF3jAGOAfEfk03XO25TLGlLxwNd0YEwC0xDqzsC2TiLwqIpVE\npArWALlFItILmGlXJmNMYOonMYwxhbHanjdi8+9U6kfrfcaY61KfugvYbHeuVN2x/lBfYGemrUAj\nY0yh1AuHd2Fd6Lf1fTLGlEr9txJwP1Yzl12ZslsvZwDdUnstVQaqYV3buzw3Xlj4CTgAJGK1y/XF\nuhCzAOuHHgWEeOIiR7pMTbA+mq3DKqZrgDZAcbtyAbVSc6wF1gMvpD5vW6YM+e7k4sVYO9+nyul+\nbhuBwXZnSpftZqwR5OuAX4GidufCurB/FAhO95zdmV7E+iO4ARgP5PeCTEux2urXYvXosuV9utp6\nidUDZwfWheNWWe1fB0wppZSPs/tirFJKKTfTQq+UUj5OC71SSvk4LfRKKeXjtNArpZSP00KvlFI+\nTgu9Ukr5OC30Sinl4/4f1UprVF1MIZAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x109371630>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t_values, t_counts = np.unique(target, return_counts=True)\n",
    "t_quantiles = np.cumsum(t_counts).astype(np.float64) / target.size\n",
    "\n",
    "plt.plot(s_values, s_quantiles, label=\"Source\")\n",
    "plt.plot(r_values, r_quantiles, label=\"Reference\")\n",
    "plt.plot(t_values, t_quantiles, '--r', lw=2, label=\"Target\")\n",
    "plt.legend(loc=5)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
